(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(*according to fastica*)
BeginPackage["W`ICA`"];
fastICA::usage="fastICA[x, minDeltaW, maxLearning]
minDeltaW:Stop criterion : Desired minimum change between \!\(\*SubscriptBox[\(W\), \(n\)]\) and \!\(\*SubscriptBox[\(W\), \(n + 1\)]\);
maxLearning:Stop criterion : Maximum number of steps allowed for convergence; m:whether return Mixing Matrix";
rowMean::usage="center row data";
k4::usage="four order cumulant(Kurtosis)";
(*Whiten::usage="";*)
(*options:*)
Epsilon::usage="Stop criterion : Desired minimum change between \!\(\*SubscriptBox[\(W\), \(n\)]\) and \!\(\*SubscriptBox[\(W\), \(n + 1\)]\)";
maxNumIterations::usage="Stop criterion : Maximum number of iteration steps allowed for convergence";
retW::usage="whether return mixing and demixing matrix, d:false";
noMean::usage="do not add mean when get ICs, d:true";
orderF::usage="function used to order ICs, orderF[Ic]=real_number, default not ordered";

Begin["`PP`"];


Options[fastICA]={Epsilon->10^-3,maxNumIterations->100,retW->False,noMean->True,(*orderF->k4*)orderF->None};


fastICA[x_,opts___]:=Module[{minDeltaW,maxLearning,ms,nbSignals,dim,xc,cov,d,e,wMat,dewMat,nx,k,done,hypTan,nw,err,sig,w,y,a1,En,rd,re,nIC,W,A,t1,t2,failureLimit,numFailures,failed,numTrys,odf,ind},

minDeltaW=Epsilon/.{opts}/.Options[fastICA];
maxLearning=maxNumIterations/.{opts}/.Options[fastICA];
ms=retW/.{opts}/.Options[fastICA];
(* output options*)
Print["use deltaW:",minDeltaW,", Maxstep:",maxLearning];

(*Print["Input Correlation Matrix= ",TableForm[Chop[Correlation[x\[Transpose]]],TableHeadings->{Table["IC "<>ToString[i],{i,Length[x]}],Table["IC "<>ToString[i],{i,Length[x]}]}]];*)

(* In fact, here ,each column is a sample,each row is a feature (as a signal) which should be independent *)
{nbSignals,dim}=Dimensions[x];

(*xc=Table[x[[i,All]]-Mean[x[[i]]],{i,nbSignals}];*)
xc=Standardize[x\[Transpose],Mean,1&]\[Transpose];

{wMat,dewMat,nIC}=Whiten[xc,Length[xc]];
(*{d,e}=Chop[Eigensystem[Chop[Covariance[xc\[Transpose]]]]];
(*See if we have to reduce the dimension*)
nIC=nbSignals-Count[d,0];
rd=d[[;;nIC]];
re=e[[;;nIC]];(* row wise*)
If[nIC!=nbSignals,
Print["Dimension adjusted to ",nIC," according to covariance matrix"]];

Print["Smallest remaining (non-zero) eigenvalue ",rd[[-1]]];
Print["Largest remaining (non-zero) eigenvalue ",rd[[1]]];
 Print["Sum of ",nbSignals-nIC," removed eigenvalues(Choped) ",Total[d[[nIC+1;;]]]];
Print[Total[rd]/Total[d]*100,"% of (non-zero) eigenvalues retained"];

(*dewhiten*)
dewMat=re\[Transpose].DiagonalMatrix[rd^(1/2)];
(*whiten matrix*)
wMat=DiagonalMatrix[rd^(-1/2)].re;*)
nx=wMat.xc;
(*each row is a feature and a signal*)

(*Print[Covariance[nx\[Transpose]]];*)
(*Print[nx.nx\[Transpose]];*)
Print["whiten OK"];

(*ICA*)
(* arg: a1,where 1 <= a1 <= 2*)(* arg: nIC*)
a1=1;
failureLimit=5;
numTrys=0;

Label[begin];
failed=False;
t1=AbsoluteTime[];
(*Random guess*)
(* nIC rows and Length[nx](num of select signal) column*)
w=Orthogonalize[Table[RandomReal[{0,1}],{nIC},{Length[nx]}],Method->"Householder"];
(*Print[w\[Transpose].w];*)
(*Print["d od w,nx:",Dimensions[w],Dimensions[nx]];*)
Print["guess OK"];

k=1;
done=False;
numFailures=0;
While[!done,
(*PrintTemporary["istep"<>ToString[k]];*)

(* using Tanh,max nongauss*)
hypTan=Tanh[a1*nx\[Transpose].w[[1]]];
(*Print[AbsoluteTime[]-t1];*)
(*nw=Flatten[(nx.hypTan-a1*w[[1]]*Sum[(1-hypTan^2)[[i]],{i,dim}])/dim];*)
nw=Flatten[(nx.hypTan-a1*w[[1]]*Total[(1-hypTan^2)])/dim];
(*Print[AbsoluteTime[]-t1];*)
nw=nw/Norm[nw];
(*Converged when nw//w *)
err={Norm[w[[1]]-nw],Norm[w[[1]]+nw]};
(*Print[AbsoluteTime[]-t1];*)
If[k==maxLearning||err[[1]]<minDeltaW||err[[2]]<minDeltaW,done=True,w[[1]]=nw;k=k+1;]
];
PrintTemporary["IC 1 computed in "<>ToString[k]<>"steps"];


For[sig=2,sig<=nIC,sig++,
(*decorrelation*)
w[[sig]]=w[[sig]]-Sum[(w[[sig]].w[[i]]*w[[i]]),{i,sig-1}];
w[[sig]]=w[[sig]]/Sqrt[w[[sig]].w[[sig]]];

k=1;
(*done=False;*)
numFailures=0;
While[True,
hypTan=Tanh[a1*nx\[Transpose].w[[sig]]];
nw=Flatten[(nx.hypTan-a1*w[[sig]]*Total[(1-hypTan^2)])/dim];
nw=nw-Sum[(nw.w[[i]]*w[[i]]),{i,sig-1}];
nw=nw/Norm[nw];
err={Norm[w[[sig]]-nw],Norm[w[[sig]]+nw]};
If[err[[1]]<minDeltaW||err[[2]]<minDeltaW,Break[],
w[[sig]]=nw;k=k+1;];

If[k==maxLearning,
(* try again? *)
numFailures++;
If[numFailures>failureLimit,
Print[sig," ",numFailures," Failed to converge, give up"];failed=True;numTrys++;Break[],
Print[numFailures," time(s) Failed to converge, try again"];
(* guess again*)
w[[sig]]=Table[RandomReal[{0,1}],{Length[nx]}];
w[[sig]]=w[[sig]]/Norm[w[[sig]]];k=1;
];
];
];(*While end*)
PrintTemporary["IC "<>ToString[sig]<>" computed in "<>ToString[k]<>" steps"];
];

t2=AbsoluteTime[]-t1;
Print["total time:",t2,"seconds"];
(*the procedure failed, should we try again?*)
If[failed,If[numTrys>failureLimit,Print["totally failed, exit"];Return[],
Print["failed ",numTrys," times, try from the beginning"];Goto[begin]],
(*ok success!*)
Print["ICs computation succeeded!"];];


(*Return["debug"];*)
(*Print["d of w,nx:",Dimensions[w],Dimensions[nx]];*)

(*mixing matrix*)
W=w.wMat;
A=dewMat.w\[Transpose];
(*Estimated IC*)(*add mean back*)
(*y=w.nx+w.Transpose[Table[(Mean[x\[Transpose]][[;;Length[nx]]]),{dim}]];*)
If[noMean/.{opts}/.Options[fastICA],
y=W.xc,
y=W.x];

(*order needed??*)
odf=orderF.{opts}/.Options[fastICA];
If[odf!= None,
(*the ICs and W,A must be adjusted at the same time*)
ind=Reverse[Ordering[Map[Abs[odf[#]]&,y]]];
y=y[[ind]];
W=W[[ind]];
A=A[[All,ind]];
];

If[!ms,
Return[y],Return[{y,W,A}]];
];


Whiten[cd_,nIC_]:=Module[{d,e,dim,rd,re},
(* if the num of signals are much larger than samples, there will be a huge cov*) 
(* choose the smaller one *)
If[Length[cd]<5 Length[cd[[1]]],
(*normal way*)
{d,e}=Chop[Eigensystem[Chop[Covariance[cd\[Transpose]]]]],
(*opt way*)
Print["using opimized method for whiten"];
{d,e}=Chop[Eigensystem[cd\[Transpose].cd/(Length[cd[[1]]]-1)]];
e=e.cd\[Transpose];
e=Map[Normalize[#]&,e]
(*Print[Dimensions[e],e,d];*)
];

(* check if lambda is 0 and find the real num of ICs*)
dim=Min[nIC,Length[d]-Count[d,0]];
rd=d[[;;dim]];
re=e[[;;dim]];(* row wise*)
If[dim<Length[cd],
Print["Dimension adjusted to ",dim," according to covariance matrix"]; Print["Sum of ",Length[cd]-dim," removed eigenvalues(Choped) ",Total[d[[dim+1;;]]]];
Print[Total[rd]/Total[d]*100,"% of (non-zero) eigenvalues retained"]];

Print["Smallest remaining (non-zero) eigenvalue ",rd[[-1]]];
Print["Largest remaining (non-zero) eigenvalue ",rd[[1]]];

(* whiten and dewhiten m*)
{DiagonalMatrix[rd^(-1/2)].re,re\[Transpose].DiagonalMatrix[rd^(1/2)],dim}
];
(*Whiten[cd_]:=Module[{},
Whiten[cd,Length[cd]]
];*)


rowMean[ss_]:=Module[{},
Standardize[ss\[Transpose],Mean,1&]\[Transpose]];

k4[d_]:=Kurtosis[d]-3;


End[];
EndPackage[];



